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It is shown in this work how the Wang-Landau algorithm can be parallehzed through the concept of 
the micromagnetic ensemble, when the Hamiltonian contains both spin interaction and the external 
field terms, and thus energy-magnetization plane is used for characterizing the density of states. 
Within this framework random walk is performed on mutually independent micromagnetic lines, 
and can thus be paralellized on a computer grid, without need for shared memory among individual 
processes. This approach pushes forward significantly the size of the systems that may be addressed 
on current computer hardware (from currently reported 42 x 42 to at least 256 x 256 in the case of 
two dimensional systems), and should turn out important for diverse studies where field dependent 
behavior is essential. 



Wang-Landau (WL) algorithm [T represents the last 
(most successful) work in a series of attempts performed 
over the last couple of decades [2-8 to generalize the 
canonical importance sampling approach of Metropolis 
et al. |9J. While importance sampling [9 traces out a 
path in the configurational space leading to equilibrium 
configurations for the particular choice of external pa- 
rameters (such as temperature and field), all of these 
novel methods are concerned with estimating the density 
of states (DOS), that is, the number g{E) of possible 
configurations available to the system at a particular en- 
ergy level E. The density of states curve is independent 
o/ temperature, it depends on the topology of the lattice 
alone, and it contains all the information necessary for 
the complete solution of the problem at hand. 

The WL algorithm is accomplished through an iter- 
ative procedure, where a random walk is performed in 
the configurational space while simultaneously augment- 
ing the density of states by a multiplicative factor / > 1, 
and incrementing a histogram of visited configurations. 
The transition probability between states with energies E 
and E' is proportional to the ratio of the corresponding 
(previously accumulated) degeneracies g{E) and g{E')^ 
and thus by construction, more probable (higher entropy) 
energy levels develop higher DOS values g{E)^ and tran- 
sition probabilities level out so that a flat histogram is 
eventually attained. The multiplicative factor / (typi- 
cally starting out with a value of In / = 1) is systemat- 
ically reduced as / ^ /4 (In/ ^ In / — 0.25) and the 
histogram is reset every time that that it fulfills a (pre- 
defined) flatness criterion (common convention is when 
80% of visited states are not lower than 80% of the his- 
togram average), until a (predefined) lower bound of / 
is reached (typically, fmin = 10~^). In short, WL is 
a highly intuitive, appealing, heuristic algorithm (with 
several somewhat arbitrary, but conventionally adopted 
parameters), which has been accepted by the scientific 
community as the state of the art, over the last decade or 
so. As at each step it (strongly) depends on the previous 



history of the ongoing simulation, WL may be regarded 
as a non-Marcovian Monte Carlo algorithm. 

If the Hamiltonian contains a field term, the density 
of states g{E^ M) is defined on the energy-magnetization 
plane, and the simulation becomes far more demanding. 
In particular, while results for zero field have been re- 
ported [1 up to sizes 256 x 256, memory requirements 
and the difficulty of convergence of the two dimensional 
random walk algorithm have up to date restricted stud- 
ies of systems in a field to sizes 42 x 42 [lOl [11]. In 
comparison, exact results for DOS functions can be ob- 
tained for the zero field case using the method proposed 
by Beale [12] with algebraic manipulation software (such 
as Mathematica or Maple) up to size 64 x 64 on cur- 
rent computer hardware, and exact DOS surfaces for the 
non-zero field case calculated using the transfer matrix 
method [13 have been reported [14 up to size 12 x 12. 
In what follows, it is shown how system size limitation of 
the Wang Landau algorithm for systems in the field may 
be significantly extended. 

Without loss of generality, let us consider here the Ising 
model with the nearest neighbor coupling J, in a uni- 
form magnetic field H (on an arbitrary lattice), with the 
Hamiltonian 

^ = -JY. S^S,-HY,S^ , (1) 

<*,J> i 

where <> denotes summation over nearest neighbor 
pairs, and 5'^ = ±1 is the spin at site i. The partition 
function of this system may be written as 

Nb N 

Z = Y.Y.9ue-''''\ (2) 

where N is the number of spins and Nb the number of 
bonds, Em = -J{Nb - 2k) - H{N - 2i) is the energy of a 
configuration having k G {0, . . . A^5} pairs of antiparallel 
spin pairs (number of "unhappy" bonds), G {0, . . . N} is 
the number of spins parallel to the field ("up" spins), and 
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Qki are the corresponding degeneracies. Setting H = 
yields microcanonical degeneracies = ^^gui while 
setting J = leads to "micromagnetic" degeneracies gi = 
Ski^ which correspond simply to the number of ways 
one can arrange i "up" spins in an N spin system, that 
is gi = (^). 

The current approach is based on the observation that 
the Wang-Landau algorithm may be performed indepen- 
dently for individual micromagnetic lines ^, as follows. 
The initial state S£{k) is prepared by flipping i spins 
starting from an ordered state with all spins up, resulting 
in k unhappy bonds, and lists of up-spin and down-spin 
indices are prepared. The initial density of states are 
set to g£{k) = 1, and the histogram values are set to 
zero h£{k) = 0. Lower index is used here to emphasize 
the fact that i is fixed, and storage space corresponds 
only to the number of bonds in the system Nb (rather 
than the product Nb x N of the number of bonds and 
the number of spins, as required by the two dimensional 
WL random walk). At each Monte Carlo step one spin 
is randomly chosen from the up-spin list and another 
from the down-spin list, and both spins are flipped to 
produce the proposal state S£{k') (therefore by construc- 
tion proposed states preserve magnetization - sampling 
is restricted within the current micromagnetic ensemble) . 
The new (proposed) state is accepted with probability 



p{si{k) si{k')) = min 



(3) 



where the density of states is modified by a multiplicative 
factor / (initially / = e) to g^ik') gi{k')^f^ histogram 
is incremented hi{k') hi{k') + 1, and the up-spin and 
down-spin lists are correspondingly updated. If the new 
state is rejected, the density of states of the current state 
and the corresponding histogram entries are updated to 
gi{k) geik) * / and h£{k) hi(k) 1^ respectively, 
while the up-spin and down-spin lists are not modified. 
When the flatness criterion is reached (e.g. of 80% of 
histogram entries being greater than the histogram av- 
erage), the histogram is reset, the density of states is 
normalized as 



i1) 



(4) 



and the multiplicative factor is reduced to / ^ / 4 . The 
simulation stops when the multiplicative factor becomes 
less than some predefined value, here fmin = 1 + A^*10~^ 
is used, since the logarithm of the density of states scales 
as N (the total number of configurations is 2^). 

In fact, it was observed that at each new level of / (af- 
ter histogram reset together with DOS renormalization) , 
the WL walker spends considerable time at the ends of 
the energy spectrum (the smallest entropy regions), and 
the histogram builds up rapidly in these regions, making 
it difficult to reach the flatness criterion. It was found 



that additional histogram resetting with density of states 
normalization at the same precision (unchanged value of 
/) is helpful in speeding up the algorithm, such that con- 
vergence with current precision is achieved in a matter 
of minutes for a single micromagnetic line for a 64 x 64 
system on a single core of a 2.4GIIz Intel 17 processor, 
and in a matter of hours for a 256 x 256 system 

The current scheme is evidently ergodic, as any con- 
figuration with i spins down can be obtained from any 
other configuration that also has i down spins, by flip- 
ping at most i spin pairs (this worst case scenario cor- 
responds to the situation when the Hamming distance 
between the two configurations assumes the maximum 
value of £). The WL random walk now becomes one 
dimensional, performed on the interaction energy levels 
corresponding to the current (chosen) micromagnetic en- 
semble. This fact, together with the reduced storage re- 
quirements, makes it possible to run much larger systems 
in comparison with the two dimensional walk. 

In Figs. [1] and [2] results are presented for L x L Ising 
model systems for L = 16,32,64,128,256 described by 
Hamiltonian ^Jot£ = N/S, 3N/S and N/2, which 
individually took up to 12 hours on a single core of an 
Intel 17 processor. 
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FIG. 1. Micromagnetic density of states i = i = iV/4, 

i — 3N/8 and i = N/2, for the LxL Ising model with periodic 
boundary conditions, for L = 16, 32, 64, 128, 256. 

The difference between the way how the systems with 
periodic and open boundary conditions approach the 
thermodynamic limit is evident on Figs, [l] and [2j and 
may be understood by considering the spin configura- 
tions that correspond to the low and high energy values. 
The low energy configurations are accomplished by clus- 
tering of up and down spins, and for a periodic system 
with Nb = 2L^ bonds two equal blocks of N/2 up and 
N/2 down spins are separated by an interface consisting 
of k = 2L frustrated bonds, such that the ground state 
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FIG. 2. Micromagnetic density of states I = I = N/A, 

i = 3iV/8 and i = N/2, for the L x L Ising model with open 
boundary conditions, for L = 16, 32, 64, 128, 256. 

corresponds to k/Ns = 1/^ (e.g. k/Ns = 0.0625 for 
L = 16). For the systems with open boundaries with 
Nb = 2L(L — 1) bonds the two blocks are separated 
by an interface of k = L frustrated bonds, so that the 
ground state corresponds to k/Ns — l/2(i^ — 1) (e.g. 
k/NB = 0.0333 for L = 16). 

On the other side of the energy spectrum both pe- 
riodic and open boundary systems with N/2 spins up 
and N/2 down spins have k = Nb frustrated bonds 
in the Neel configurations, so that k/NB = 1. For 
lower values of £ the high energy levels correspond 
to scattering of up spins, such that none are neigh- 
bors of each other, yielding k = M frustrated bonds, 
and therefore the upper energy bound is attained at 
k/NB = U/Nb {k/NB = N/2Nb,N/Nb,3N/2Nb for 
£ = 7V/8, 7V/4, 3Ar/8, respectively). For periodic bound- 
aries N/Nb = 1/2 and the density of states curves on 
Fig.[l]end at 0.25, 0.5, 0.75 for £ = N/8, N/4, 37V/8, while 
for open boundaries N/Nb = L/2{L — 1), and the lim- 
iting upper energy bound levels are being gradually ap- 
proached as k/NB = L/4(L - 1), L/2(L - 1), 3L/4(L - 1) 
for £ = N/8, N/4, 3N/8, respectively. 

Finite size scaling of the DOS curves evidently re- 
quires adjusting both axes (the dimensionless energy U 
and the DOS function Ski = ^^9k£{L) scales), and does 
not appear to be a straightforward matter (studies in 
this direction are under way, and any conclusive results 
shall be reported elsewhere). On the other hand, gen- 
eral scaling behavior may be inferred by observing only 
the DOS function values Sc{L) at the center of the en- 
ergy magnetization plane {k/NB = 0.5, £/N = 0.5), 
as a function of system size. In particular, the quan- 
tity AS{L) = ln2 - Sc{L) versus ln(l/L) (where In 2 is 
the limiting maximum entropy value) demonstrates lin- 
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FIG. 3. Difference of the scaled entropy value and the 
maximum value In 2 at the center of the energy magnetiza- 
tion plane, as a function of inverse linear size, for the L x L 
Ising model with open and periodic boundary conditions, for 
L = 16,32,64,128,256. 

ear behavior for both periodic and open boundary con- 
ditions, as shown in Fig. [3] Linear behavior of AS{L) 
versus ln(l/L) observed in Fig. [3] implies scaling form 

Sc{L) = ln2-^, (5) 

where values obtained by regression are a = 1.252 and 
P = 1.778 for open boundaries, while for periodic bound- 
ary conditions parameter values a = 1.088 and /3 = 1.758 
are obtained. For other points on the energy magnetiza- 
tion plane one may expect similar scaling behavior 

Ski{L) = Ski{oo) - (6) 

where Ski{oo) represents the corresponding entropy value 
in the thermodynamic limit. 

It is seen that the micromagnetic DOS functions can 
be determined independently of each other using the 
current implementation of the WL algorithm, for very 
large systems. Nevertheless, determining the complete 
set of micromagnetic lines that comprise the full DOS 
surface above the energy magnetization plane remains 
a formidable problem in terms of computer resources re- 
quirements, with a (rough) estimate of 200000 core-hours 
for the total of 32768 micromagnetic lines of a 256 x 256 
system (for models with up down symmetry, or double 
that number for asymmetric models), at the current pre- 
cision of / = 1 + TV * 10~^. On the other hand, as cal- 
culations for individual micromagnetic lines are fully in- 
dependent of each other (no parameter region border ad- 
justment, or relative normalization is needed), they may 
be performed on a geographically distributed computing 
grid, and the task may be regarded as rather demanding, 
but feasible. 
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FIG. 4. Density of states for the 32 x 32 and 64 x 64 Ising 
model systems with periodic boundary conditions, calculated 
using the current implementation of the WL algorithm. 

The full DOS surfaces for L x L systems with periodic 
boundary conditions, for L = 32 (obtained in a matter 
of hours on a single Intel 17 processor) and for L = 64 
(obtained in a matter of days) are shown in Fig. [4} In 
accordance with results displayed in Fig. [l] it is seen 
on Fig. |4] that the DOS surfaces are rather similar, the 
largest differences being observed in the region of low 
energy {k/Ns ~ 0) and low absolute magnetization val- 
ues {i/N ^ 0.5). The depression of the DOS surface 
in this region, surrounded by symmetric ridges that join 
smoothly with energy increase, is the signature of the 
second order transition. 




FIG. 5. The exact density of states Sk = Ingk/Ns for the 
64 X 64 Ising model with periodic boundary conditions, to- 
gether with the results obtained by summing (over i) the 
micromagnetic curves obtained through application of the 
Wang-Landau algorithm. On the right side the difference be- 
tween the exact and Wang-Landau algorithm results is shown. 

While the exact solution for the Ising model in the 
field is not known and there is no known exact result 
with which the DOS surfaces of Fig. [4] can be compared, 
one further test of validity can be made by summing the 
density of states over micromagnetic variable for each 
energy level, and than comparing the result with the ex- 
act density of states in zero field, obtained through the 
method proposed by Beale [12] using algebraic manipula- 
tion software (such as Mathematica or Maple). Results of 
this comparison are displayed in Fig. |5j where no visible 
difference is seen on the scale of the graph. The difference 
between the exact dos functions and the WL estimate 
shown on the right hand side of Fig. [5] is seen to be well 



below the value ^/\nf = \^N10^ - 0.002 [15] (except at 
the very ends of the energy spectrum range), which may 
be attributed [15] to multiple simulations performed for 
distinct micromagnetic ensembles (different values of ^), 
that were used here to compose the microcanonical en- 
tropy curve. 

In summary, in this work it is shown how the Wang 
Landau algorithm may be optimized for systems with 
a field dependent Hamiltonian, by running independent 
micromagnetic ensemble runs. Updates are performed by 
conserving magnetization (number of spins parallel to the 
field), and this reduction of the dimension of the param- 
eter space brings about considerable advantages. First, 
memory requirements are reduced from Nb x TV (where 
Nb is the number of bonds, and N number of spins in 
the system) to Nb- Second, the micromagnetic ensemble 
runs are fully independent of each other and may be per- 
formed in parallel on a geographically distributed grid. 
And finally, convergence of the algorithm is observed on 
much larger systems as compared with a two dimensional 
WL random walk (pushing the limit from 42 x 42 reported 
field results [10l[TT] to at least 256 x 256). This approach 
is independent of the details of geometry and interac- 
tions (although examples presented here were performed 
on a nearest neighbor Ising ferromagnet), as long as the 
uniform field term is present in the Hamiltonian. 
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